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We study the limitations for entanglement due to collisional decoherence in a Bose-Einstein 
condensate. Specifically we consider relative number squeezing between photons and atoms cou- 
pled out from a homogeneous condensate. We study the decay of excited quasiparticle modes 
due to collisions, in condensates of atoms with one or two internal degrees of freedom. The time 
evolution of these modes is determined in the linear response approximation to the deviation 
from equilibrium. We use Heisenberg-Langevin equations to derive equations of motion for the 
densities and higher correlation functions which determine the squeezing. In this way we can 
show that decoherence due to quasiparticle interactions imposes an important limit on the degree 
of number squeezing which may be achieved. Our results are also relevant for the determination 
of decoherence times in other experiments based on entanglement, e.g. the slowing and stopping 
of light in condensed atomic gases using dark states. 

PACS numbers: 03.75.Fi, 03.67.-a, 42.50.-p 
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I. INTRODUCTION 



Squeezing and entanglement in coherent many particle systems have become important issues in experiments 
with Bose-Einstein condensates (BECs) as well as in the theoretical work related to them. Achieving squeezed 
■ and/or entangled states of a large number of ultracold atoms would have an enormous impact on the research field 
which includes quantum information technology (QIT) [Q, the development of new time standards j^j and the test 
of quantum state reduction ||. In particular, whereas techniques for entangling photons Q are already available 
for the transmission of quantum information ||, schemes for the local storage and processing are still missing. A 
£Sj ■ promising step in the desired direction represent the experiments which allow for photon pulses being slowed and 
even stored for a considerably long time in cold and ultracold, condensed atomic gases || (for a recent review 
cf. These experiments are based on electromagnetically induced transparancy [|| and use the possibility of 
l dark states, which for the case of many-particle systems involving a condensate may be described as entangled 
^Ij ■ states of atoms and photons. 

In Q we have studied the limitations of relative number squeezing in entangled states of atoms and photons 
produced by Bragg-scattcring of light off a condensate. Entanglement between photons and atoms has been 
qh proposed and studied by different authors before . We took into account interactions of the condensate atoms 
with non-condensed atoms to the degree that the quasiparticles resulting from the diagonalization of the quadratic 
part of the Hamiltonian remained free. We have shown that the amount of squeezing in such systems underlies 
important limits imposed by the collisional interactions of the atoms. In this paper we would like to extend this 
work and study the effect of quasiparticle interactions which lead to collisional damping of excitations. We will 
find additional constraints on the squeezing considered in |j| . 

Moreover, we will consider the case where the non-equilibrium matter wave excitations consist of atoms which 
are in a different internal state than the condensate atoms. This is e.g. relevant in the case of a dark state involving 
a condensate, a classical coupling laser and quantized probe and atomic modes, which allows to slow and stop 
light pulses in atomic matter JlT| . The quantized atomic mode is in general in a different internal state than the 
classical (condensate) mode and the decoherence time is governed by the spin-conserving collisions between the 
atoms of the different modes. In this paper we derive the collisional damping rate of an excited atomic mode for 
such a system with two internal degrees of freedom, in an analogous way as in the single level case. In pd| ] we have 
studied the limitations due to this collisional damping of the light delay and storage times in experiments based 
on electromagnetically induced transparency. 

Mean field theory and its extensions have been widely used to describe stationary properties of homogeneous and 
trapped dilute, weakly interacting Bose-Einstein condensates at finite temperatures |l^-|l4|]. In these approaches, 
the non-condensate component is treated as a static thermal bath and its dynamical coupling to the condensate 
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oscillations is neglected. This implies that damping of out-of-equilibrium-excitations is not accounted for. Time 
dependent mean field approaches have been introduced Jl5|— |l8|] to describe the damping in trapped condensates 
at finite temperatures Jul . The calculated damping rates agree with the results from perturbation theory jl9],|o| . 
In [^l],^2| damping rates in trapped condensates at finite temperature have been calculated explicitly. In order 
to describe the dynamics of a decaying excitation consistently at large times it is necessary to take into account 
the connection between dissipation and fluctuations. To this end it is convenient to use the Heisenberg-Langevin 
approach (cf. e.g. p3||) to treat the time development of the system p4] , p2| . In this paper we would like to recall 
this procedure as applied to the dissipation dynamics of condensates with atoms in a single internal state and 
extend it to the treatment of the two internal levels case. Moreover, the Heisenberg-Langevin equations will enable 
us to derive the time evolution of expectation values of products of two and more quasiparticle operators which 
we need to study the squeezing described above. 

In the basic, quadratic approximation to the Hamiltonian of the system, including a quartic collisional contact 
interaction, mean field theory yields the ground state wave function and the spectrum of elementary excitations, 
which may be understood as non- interacting quasiparticles. In this theory, the Heisenberg equations of motion of 
the Fock operators assume the form 

kit) = -mPi(t), (i) 

where Pi{t) is the annihilation operator of the quasiparticle mode labelled by i, which satisfies Bose equal-time 
commutation relations with its hermitian conjugate, 0i(t), $j{t)) — Sij. uii is the corresponding eigenfrequency. 

Beyond this free field theory the atomic collisions induce interactions of the quasiparticles. The purpose of this 
paper is to consider the impact of these interactions on the dynamic behaviour of a BEC system which is initially 
in a state out of thermal equilibrium. Specifically we study the decay behaviour of a single quasiparticle mode 
whose initial population is different from the equilibrium one. As this quasiparticle mode is interacting with the 
bath of other modes the system will evolve back to a stationary equilibrium state. The dynamics of such processes 
may be described using a time-dependent mean field approach based on the finite temperature extension of the 
Gross-Pitaevski equation jl6],[L8) . This allows to describe oscillations and damping of the condensate mean field 
which initially is assumed to be slightly displaced from its stationary form. 

As an alternative to the mean-field approach one may adhere for the first to operators and find the equations of 
motion for the quasiparticle operators. In the basis of the non-interacting quasiparticles of the quadratic theory 
these will be in general coupled non-linear equations. For the physical situations we are considering here, we are 
interested in approximations which allow to write the equation for the excited mode q in the form 

$ q (t) = -i(u q + 6w q -^ q )$ q {t), (2) 

where 5uJn is an energy shift and 7 9 a decay width. As was shown using the extended time dependent mean-field 
theory the decay width arises from the coupling of the excited mode to the bath of the other modes 

which are assumed to be in thermal equilibrium and to have a broad range of suffiently closely spaced frequencies. 
Eqn. (|^) may be used e.g. to calculate the time dependence of the mean population of mode q, which results as an 
exponential decay with rate 7 9 . It is also clear, however, that the operator which evolves in time as described by 
(||) does not obey the equal time commutation relation, but commutes with /3j(i) for t — > oo. Moreover, according 
to (||) the population of the mode q decays to zero, i.e. not to the correct equilibrium value at finite temperatures. 
These problems are absent in the mean-field approach, where only the deviation from the stationary condensate 
field decays and where no operator ordering problems occur. 

It is well known in quantum optics that these problems are resolved by the Heisenberg-Langevin theory of 
dissipative dynamics. Let us briefly summarize its results. The equation of motion (||) for /3 q contains a further 
term. In the Markov approximation, the equation may be written in the form 

4(t) = -*K + <H - hi) Kit) - F q (t), (3) 

where F q is a Langevin quantum noise operator which has a zero mean value. The Langevin operator in the 
Markov approximation obeys the commutation relation 

>,(t),#(f)l = lq Sit-t'). (4) 
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This ensures that the equal time commutation relation for (3 q {t) is conserved in time, 



q (t)Jl(t)} = 1. (5) 

In the Markov approximation the noise operator is delta-correlated in time, i.e. the two-time correlation function 
is given by 

{F^t)F q (t'))= 7q n q S(t-t% (6) 

where n® is the equilibrium value of the mode population. Eqn. (|J) means, that the random forces described by 
F q have a short memory. The equation of motion for the quasiparticle number operator n q — (3^ q (3 q becomes 

h q {t) = - lq [n q {t)-n Q q ]-F' q {t), (7) 

where F' q is the quantum noise operator defined by 

F' q = p\P q + Ffp q -0\F q + Pfp q ) , (8) 



n° 



which by definition has a zero mean value. To prove eqn. (jjp, it has to be shown expicitly that (f3^P q + P^$, 
Eqn. (R) shows that the quasiparticle population of mode q decays, as expected, exponentially to the equilibrium 
value rr q . Before the population can decay to zero the fluctuations introduced by the coupling to the bath modes 
force it to settle at the thermal equilibrium value. 

Our paper is organized as follows: In section |n] we review the Heisenberg-Langevin approach to quasiparticle 
damping for the single internal level case. We use a linear response treatment of the operator time evolution 
analogous to the one used in Jlj|[l8| for the mean field fluctuation and derive both (Q) and (^) in the Markov 
approximation. In section |lj we extend the approach to the case of two internal levels and derive the Beliaev- 



and Landau-type decay widths as well as the corresponding noise operators. In section [V we study time evolution 
equations of the form (Q) for normal ordered products of two and more quasiparticle operators. We introduce 
a quantum regression theorem which allows to use the basic equation (^) for the operator to calculate the time 
evolution of the operator products. We use this to calculate the implications of quasiparticle damping for the time 
evolution of the atom-photon squeezing produced by Bragg scattering of photons off condensates as described in 
||. Section |v| contains our conclusions. 



II. HEISENBERG-LANGEVIN EQUATION FOR SYSTEMS OF ATOMS IN A SINGLE INTERNAL 

STATE 

The Hamiltonian of the interacting system of non-relativistic bosonic atoms in the presence of an external 
trapping potential V(r), in the grand canonical description, has the form 



K = H-^N = J d\ j^(r,i) r-^!+y(r)-^(r,t) + |^( r , t )^t( r>t )^( r>t )^( r>t )|. 



We have chosen the contact potential approximation for the interatomic potential, U(r) — gS(r), where g is given 
in terms of the s-wave scattering length by g = Ai:h 2 a/m. The field operators ip(r,t) and $t(r,i) obey the usual 
equal-time commutation relations. Using these, the equation of motion for the field operator follows to be 

ih^(r,t)= + ^(r,t) + ffVi t (r,t)Vi(r,t)Vi(r,t). (10) 

In the usual treatment of the time evolution of a Bose condensate away from, but close to equilibrium a time 
dependent condensate wave function is defined by 

$(r,i) =<$(!■,*)), (11) 

where the avarage ( • • • ) is with respect to the non-equilibrium initial state of the system. The field operator is 
then decomposed into a condensate and a non-condensate component, 
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^(r,t) = *(r,i)+$(r,*), 



(12) 



and the non-condensate part is defined by (t/>(r, t)) = Eq]. The condensate part is usually taken to be a 
c-number, but we will keep its operator character as far as its non-equilibrium contribution is concerned: 



®(r,t) = $ o (r) + 0(r,t). 



(13) 



Here $o(r) = (' ( M r ))o is the stationary expectation value w.r.t. the equilibrium state of the BEC and assumed 
to be real. In the usual mean field treatment |l6| |l8| the condensate part $ is replaced by its non-equilibrium 
expectation value <E> = (<&), and therefore <j> by the fluctuation field 5$ = (<f>), which is assumed to be small 
compared to $o. In our further procedure we will also assume that the non-equilibrium expectation value of 4> 
remains small compared to $o, in order to justify an analogous linear response approach to the time-evolution of 
the field operator </>. 

We now define the fluctuation operators of the normal and anomalous densities as 



7,0 



5h = tp'tp - 
Sm = 'tfnp — to , 



(14) 
(15) 



where we have, as we will do in the following, suppressed the space and time arguments of the operators and 
densities. The non-equilibrium expectation values of the operators ( |l4l , [l5|) are the fluctuation functions Sn, dm 
defined in p6[ , which describe the deviation from the equilibrium expectation values n° = (ip^'ip)o, to = (^!^>)q. 
Inserting eqns. into © we obtain, to linear order in </>, 



if i 



dt 



(i> + = (k + 2g [n + n°] ) $ + g [n Q + m°] ft 



+ (A' ( J 



no + n + $ ^ T + $0-0 



n + m° + 2$oV' 



2dfi + Sm 



(16) 



where we have defined Ko = — V 2 /(2to) + V — fj,, and where no = |*o| 2 is the condensate density. In deriving 
( |l6| ) we have also used that the equilibrium condensate wave function <E>o satisfies the generalized stationary 
Gross-Pitaevskii equation 

(17) 



(18) 



(K +g[n + 2n° + m ]) $ = 0, 
and we have neglected the fluctuation operators 

Sffiijrip) = ftW> ~ 2"°^ - m°ft. 



In the case, where the system is in equilibrium, the r.h.s. Jl6| ) reduces to its first line. Then, eqn. ( [l6| ) may be 
solved by introducing the Bogoliubov expansion 



0(r, t)=£[«i(r) h (t)+v*(r)pl(t) 



(19) 



of the field operator in terms of quasiparticle mode operators. In equilibrium, (16) is linear in the field operators, 
one has fli(t) = /3j(io) exp(— iu>i(t — to)) and the mode functions u,, Vi are the solutions of the Bogoliubov-de 
Gennes equations pfil] 



£(r) M(r) 
M*{r) -£*(r) 



Ui(r) 
Vi(r) 



= hu>j 



Ui(r) 
Vi(r) 



(20) 



where £(r) = K (r) + 2g[n (r) + n°(r)], M(r) = g[n (r) 
by 



(r)]. A normalization of the mode functions Ui, Vi 



d\ [uJ(r) Uj (r) - v*{r)v 3 {r)} = 6 K 



(21) 
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ensures that the quasiparticle operators /3,-, flj satisfy bosonic commutation relations 



Using ( |19[ ) we can expand the fluctuation operators (jf 4j), ( |15| ) in terms of the operators 
(the notation using / and g follows jl6| ), and find 



(22) 
(23) 



(24) 
(25) 



We now proceed with solving eqn. (|Tq ) in the non-equilibrium case. For simplicity, we assume that the initial 
state is, in the basis of the quasiparticle modes defined by (|l9|), a quasiparticle number state. We assume that 
a single quasiparticle mode is populated much stronger than the thermal equilibrium population, n q (t = 0) = 
<$(Q)&(0)> » n° = = [exp(^/fc s T) - l]- 1 (< ^ «). 

We write the fluctuation operator in the form 



(r,t) = u,(r)&(t) + «:(r)#(t). 



(26) 



In order to justify the linear response treatment, the particle population is still assumed to be much lower than 
that of the condensate mode, ( <^(r, Q)0(r, 0) ) = J d 3 r ([\u q (r)\ 2 + \v q (r)\ 2 ]n q (t) + \v q (r)\ 2 ) < n . 

Note that in (^6|) we assume $ to be defined by ([l9), but with the sum over i excluding q. Following (|l4|) and 
(H) the sums over i and j in the expansions ( p4[) and J25|) then also exclude g. However, we will include the terms 
proportional to tpcj), ip^cj) an d in (|T^) into the term g<&Q[25n + Srh] which means that in (p4{), ( |25| ) the sums 
over i and j also include g, except the combination i = j = q, which is neglected in linear response. 

Now, inserting (|l^) and ( ^6|) in (|I^) and using ( ^ii| ) and d|l]), we find the equation of motion for (3 q (t): 



(27) 



where 



4 • • — 



5, 



and 



-4, 



y d 3 r$ + v q B i:j ] . 

= u*u 3 +V* Vj +V*Uj, 



Bij = muj + mVj + ViUj, 
dj = ViVj + mvj + ViUj. 



(28) 



(29) 



The time evolution equations for the density and pair operators /y , gij may be derived in an analogous way as in 
@. Inserting @ and @ into (|), using @, @, (jT|), (]l9|), and @, and keeping only terms up to linear 
order in the fluctuation operators, one finds: 



+ g j d\ $ 20^^ + ftipip + h.< 
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+ g I d\ 



4 5n fttp + dm ft ft + dm* 'ijj'tp 
5(ftfti>ft + $ (s(ftftft + 5(ft$ft 
(2n°+m°*)0 + h.c 



where 



5(ftft^ft = ftft-^-ip — 4( ) — ( r^^t ) ^ _ ( ^ ) ^t^t^ 
5{ftftft) = ftft$ - 2ti°ft - rh°*4>, 
5n = (ftip)-h = (5h), 
8m = ( ipip ) — m° — ( (5m ) . 



(30) 



(31) 
(32) 
(33) 
(34) 



Then, using (pL9|) , (|22|), ( |23j ) the equations of motion for /y and read (i, j 7^ g) 



at 



at 



5i3 



-j(wj+w i )5y--2^ S r(l + n?+n^) / d 3 r$ 



4i« 



(35) 



(36) 



In deriving the time evolution ( p5| , p6[ ) of the fluctuation operators from (|3(j), the terms proportional 8fi and dm 
which describe the coupling to the fluctuations of the bath modes have been neglected. Also the terms arising 
from the fluctuations 5(ftfttp t ft, 5{ftftft have not been included. Note also, that in the cases i = q or j = q, 
for which the differential equations (|3^) and ([36]) take the same form, further terms in K have to be taken into 
account. These terms are of second order in <j>, ft and have been neglected in (|30|). We will in the following assume 
these special cases to be included. Equations (|35| ) and ([36]) are then the analogues of eqns. (26) and (27) of |]l6| . 
We can integrate (|35| ) and (|3^) formally and insert the solutions back into (|27j). 

12 e i(wi-w s )(t-t') 



-2(g/hfJ2 I dt{2(n°-n°)\A qA 
H J 

+ (1+ n? + n°) [|B 9 ,« I 2 e-*^^^*-*') - |5 9)ii | 2 e *(««+"i)(*-*')l J, 



(37) 



The terms involving the time integrals then lead, in the Markov approximation, to an energy shift 5u> q and a decay 



width j q . The energy shift may be calculated in the Popov approximation (7 



0) 



where it vanishes for 



<j — > 0, i.e. leaves the spectrum gapless as required by the Hugenholtz-Pines theorem [E8UT3]. For temperatures 
kgT <C /i, however, the anomalous density m° can not be neglected in the calculations. A self consistent solution 
of the Hartree-Fock-Bogoliubov equations including rh , together with a renormalization of the coupling g results 
in a gapless excitation spectrum |29]| . 

We are here primarily interested in the damping of the excited quasiparticle mode. The decay width 7 g follows, 
analogously to Jl6|] to be 

i q = iL{q) +iB{q), (38) 

7l and 75 are known as the Landau and Beliaev decay widths M,naE0|3G] and are found from (p37h to be: 
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IB 



(?) = 4^(.g/?i) 2 J] |S 9)ii | 2 (l +n? + nj) S(u q - Wi - u>j). 



(39) 
(40) 



To derive ( j39| , |40| ) we have chosen the Markov approximation. This means that we assume that the sum over i and 
j over the appearing products of coupling functions (p8|), multiplied by the energy exponentials, give functions 
of r = t — t' which are sharply peaked around r = 0. The conditions for this assumption to hold are generally 
that the energies of the bath modes are closely spaced compared to the energy of the decaying mode and that 
the couplings are sufficiently well behaved at large energies. The time dependence of the sharply peaked function 
may then be approximated as Markovian, i.e. by a delta-distribution, and integration over time yields the energy 
conservation delta-distributions in ( |39| , ^o| ). 

The first term on the r.h.s. of ([37|) is the Langevin noise operator corresponding to the decay of the quasiparticle 
mode. Hence, the equation of motion for (3 q may be written in the form (||), where the decay width j q is given by 
(|3^)-([io|) and the Langevin noise operator is 



/y(o) 



+-)%(0)} 



(41) 



Using this expression we will now show that in the Markov approximation F q obeys the commutation relation 
and its two-time correlation function is given by (|^). 
To this end we need the commutation relations 



&i(0), $1,(0)1 = (6uS jk + faSj,) (1 + n° + n°) 



(42) 
(43) 



In deriving these we have neglected on the right hand sides terms proportional to the fluctuation operators fy. 
Using (E2IE3]) and the fact, that mixed commutators of and gti only yield fluctation operators g^ we arrive at 



F q {t)e^\Fl{t')e-^'] = 2( 5 A) 2 £ {2(n? - n°) \A qAj \ 2 e 



i{uj q -\-uji—ujj)(t — t ) 



+ {l + n° + n°) \\B q ^\ 2 e^-^-^-^ 



\B„ 



|2 i(uj q +u>i+Uj)(t— {') 



(44) 



In the Markov approximation the sum over i,j is replaced by a function proportional to S(t — t'). Integration of 
flU) and comparison with (p8|)-(^0|) then shows that the coefficient of this delta-distribution is given by 7 g , as 
stated in (||). 

Eqn. (|^) may be proven in a similar way, using 



fl j (0)fki(0)) = S ik 5 jl n j (n° i +l), 



(45) 
(46) 



In deriving these relations we have made use of Wick's theorem. Taking into account still the respective energy 
matchings in the different terms, the population factors become 



5{uj q + LJi - a,) n°(n° + 1) = S(u q + uJi - to,) n°(n° - n% 

8{w q - UJi - ujj) n°n° = 8(ui q - Wj - ojj) n° q (l + n° + n°). 



(47) 
(48) 
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where we used the Planck distribution n\ — {(3j(3i)o — [exp(hui / IsbT) — 1] 1 . Together with ( ^3| ) we arrive at 

(F^{t)F q {t'))e-^ t - t "> = n° q \P q (t),F^(t')\ e -W*-f ). (49) 

where the r.h.s. is given by (|44|). This yields, in the Markov approximation, the identity (0) for the correlation 
function. 

III. HEISENBERG-L ANGEVIN EQUATION FOR SYSTEMS OF ATOMS IN TWO INTERNAL 

STATES 

In the following we would like to extend the procedure reviewed in sec. [H] to the case with two internal atomic 
degrees of freedom, where the atoms in the excited mode q are in a different internal state than the condensate 
atoms. Such situations are important e.g. in the context of the dynamics of dark states involving Bosc-Einstein 
condensates which are the basis for experiments using electromagnetically induced transparency |8| for the slowing 
and stopping of light pulses in atomic matter Q . The decay of the dark state induced by atomic collisions limits 
the maximum time for which a light pulse may be stored in collective atomic excitations of condensates JlTJ . 

The interactions which cause the quasiparticle decay will be mainly collisions between atoms in the excited mode 
and condensate atoms, such that they do not change the respective internal states. Let us call the two internal 
levels of the atoms b and c, the Bose-condensed atoms being in internal state | b ) . The Hamiltonian in the grand 
canonical ensemble then has the form 

K = H-fiN = K b + K c + K bc , (50) 

Kb = J d\ + 9 -f WMb) , (51) 

K c = J d\ (k q , c + LJ cb ) tp c + 9 f » (52) 

K bc = J d\ g bc ^l^\^ c ^ b , (53) 

where Kq^ v — — h 2 V 2 / (2m) + V v (r) — fi, v = b,c, and ui cb = uj c — uj b is the frequency difference of the internal levels. 
t^,„, /i, v — b,c are the collisonal coupling constants for the respective processes being related to the respective 
scattering lengths by g^ v = 47r?i 2 a ;j „ jm. As for the single level case we split the field operators for the atoms in 
level b into mean-field and fluctuation parts, 

&(r,t) = $ (r)+&(r,t), (54) 

and the operator for the atoms in c into fluctuation and excited modes: 

i> c (r,t)=$ c (r,t) + $ c {r,t), (55) 

With these, we find the equation of motion for the excitation operator <p c in an analogous way as in the single 
level case: 



«?jj^c = (-Kq,c + 2# cc n° + g bc [n + hi] ) <j> c + g cc m° c (j>. 



5h bc + Sm bc 



(56) 



where ttP v = (V'J' ! / ; y)oj m ° v — (4'u4'u)o are the mean equilibrium density and anomalous density of the non- 
condensed fraction in v = b,c. The fluctuation operators 6n bc = 'tpl'ipc, Srh bc — ip^c are defined as in ( |l4| , |l5| ), with 
their equilibrium mean values though being zero. 

In order to derive from (pq) the equation of motion for the Fock operator of the excited mode we use the 



decomposition (19) for tp b and a particle mode decomposition for atoms in the internal state c: 

Mr,t) = J2$c,i(rMt)- (57) 
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Here the mode functions tp Cy i are solutions of the equation 

£ c (r)V^(r) =hu?i> c ,i(r), (58) 

where C c (r) = Ko iC (r)+u}cb+2g cc n^,+gbc [no + n^] , and where we have chosen the approximation M c {r) = g cc fhP c = 
0. Note that the eigenfunctions then are not quasiparticle modes since the zero mode in c is not macroscopically 
populated. The equation of motion for c q then reads 



Kit) 



(59) 



where 



A bc — 



B 



be 



(fir $ 



4>c,q( U i + V i)tit,j 



and 



fij ~ A &j > 



PiCj 



~bc 



(60) 
(61) 



(62) 
(63) 



In order to solve ( p9[ ) we need, as in the single level case, the time evolution of the operators fij,g b j- To this end 
we write down the relevant terms of the Hamiltonian K which contribute to the commutator with these operators. 
Using (ph and pi) as well as f>7|) and (p>|) they read 



K c + k bc = £ ftwic\ci + g cc J d\ 8$lft c ipjp c ), 

(fir 3>o i>t$t^ c + 4>\tictib + h.c. 



+ gbc 
+ g^ 
+ g^ 



(fir 



d'r 



Sh b iplipc + Sh c 'tpl'tpb + (Sh bc fyipc + 8fh* bc dctib + h.c.^) 



+ h.c. 



where 



S(^Jl^ b ) = ftJl^b - ( ft c ij>c ) tiltib - ( tiltib ) ftjc - (( tittib > tittic - ( tiW b ) tictib + 

5(4>i4>i'4>b) = ipti>ti>b - n° c tib, 



h.c 



Sn c = (til^. 



cr'C 

Itic), 



(64) 



(65) 

(66) 
(67) 



and 5h bc = (Sh bc ), 5m bc — {8fii bc ). The operator 5(tpltplip c ip c ) and the functions 5fi b are defined as in (j3lj) and 
(|33|) respectively. From the terms on the right hand side of (p3) only the quadratic one and the ones in the second 



line will contribute to the leading order time evolution of /y and g^. The equations of motion for these operators 
are then found to be 



-g$ = -i(u>3 + u»)g% - -g bc {nl - n*)B b q %c q {t), 



(68) 
(69) 



where we have neglected terms of higher order in the fluctuation opertors and used <f> q = ^ c , q c q . n\ and rij are the 
equilibrium populations of modes i and j in internal states b and c, respectively. Analogously to the single level 
case we then obtain, in the Markov approximation, the equation of motion for c q in the form (|^). We find that 
the decay rate 7 9 = 7l(<?) + Jb(q) has Landau- and Beliaev-like contributions (138) which read 
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7i (g) = 27r(W*) 2 ]T |^y| a (n? - n°) + - (70) 
7B (g) = 27r( 56c /^) 2 J] |B^| 2 (1 + n\ + S(cj q - Ui - cjj), (71) 



i j 

The Langevin noise operator is given by 



F q {t) = £<&c£ {4$ e-<(W? " wi) */«?(0) 

Sh*e- i K+"*)*g|j : (0)}. (72) 



9,y 



The proofs that (72) obeys the commutator (^) and eqn. (||) for the two-time correlation function are conducted 
in analogy to the single mode case, eqns. (^)-(^TJ). 

Our result ( |70|j7l|) for the damping rate j q — Jl(q) + 7b(<?) may e.g. be used to calculate the loss of coherence in 
a dark state used to delay and store light pulses. In jy] we have calculated the damping rate for a homogeneous 
condensate and found that the Beliaev rate for small momenta q of the decaying mode is proportional to q 5 as in 
the case where the internal atomic state is the same as that of the condensate atoms. The rate differs only by a 
constant numerical factor: 7c(<z)t=o ~ hq 5 / 96mirno compared to 7c(?)t=o ~ 3hq 5 /320mirno in the single level 
case (assuming abc — abb)- 



IV. DISSIPATIVE DYNAMICS OF HIGHER CORRELATION FUNCTIONS: ATOM-PHOTON 

SQUEEZING 

We would now like to apply our formalism in order to calculate the dissipative dynamics of expectation values 
of two and four quasiparticle operators. This will enable us to compute the time evolution of the relative number 
squeezing between atoms and photons produced in photon scattering off condensates. Such schemes have been 
proposed in [jL0[ . In extending this work in JTJ] we have considered the production of atom-photon squeezing using 
a finite homogeneous Bose-Einstein condensate of interacting atoms at zero temperature, illuminated by laser light 
with given frequency f2j detuned slightly from an atomic transition, cf. fig.[j]. 




FIG. 1. Bragg transition between different modes of the center of mass motion of atoms in the internal ground state. The motional 
structure of the excited state may be neglected. 



The laser photons with momentum ki predominantly interact with the condensate atoms (k = 0) and may 
transfer them in a centre of mass mode with momentum Ak (see fig. 0) |3^] . Such processes yield a state with 
relative number squeezing in a particular mode pair, i.e. the variance of the relative occupation number is suppressed 
as compared to the classical case. We computed the squeezing parameter 

& := [A(n a -n b )} 2 /(h a + fib}, (73) 

which measures the relative number variance [A(n — rib)] 2 — ( (^o ~ fib) 2 ) ~ ( n a — fib ) 2 in relation to its classical 
limit given by the sum of mean occupations (fi a + hb) = ( a^a+b^b ) , of the scattered photon mode and the recoiling 
atom mode. (The index '3' of £3 refers to the analogy of the relative number squeezing with spin squeezing. £3 
measures the variance of the 3-componcnt J3 = h(h a — fib) compared to half the total value J/2 = \{n a + hb) of 
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an abstract angular momentum. This is defined by the Fock operators a, 6, using the Schwinger representation of 
angular momentum. For the issue of spin squeezing and the definition ( |73| ) cf . |3l| . It will be instructive to study 
also the squeezing of J\ = ^{a)b + b'a) and J 2 = ^{d^b — b'a) for which we define := [A Ji] 2 /{ J/2), i = 1, 2. 

In the classical (uncorrelated) limit £3 = 1, whereas for maximum squeezing we have £3 = 0. There are 
two reasons why such perfect squeezing is spoilt in an experiment: Collisional atomic interactions and photon 
(rescattering) processes different from those which produce the squeezed state. In [51 we focused on the limitations 
to the maximum squeezing due to those atomic interactions which are taken into account in the basic Hartree- 
Fock-Bogoliubov approximation, in the Hamiltonian quadratic in the particle operators. There we have shown 
that "unwanted" collisions amongst the particles present profound limitations to the achievable squeezing, with 
correlations between atomic modes in the initial state being the main limiting factor at small times. 

Here we would like to extend our work to include interactions between quasiparticles described by the third order 
terms in the Hamiltonian considered in section |n|. The main effect at this order of the approximation is collisional 
loss from the excited modes which will induce decay of the population of the recoiling atomic mode. This will also 
affect the time evolution of the correlation functions (An a ) 2 , (An^) 2 , and ((n a hb)) = {h a hb) — (fi a )(hb), which 



enter the squeezing parameter (73). We will neglect energy shifts at this order. 



In order to be able to compute the dynamics of the occupation numbers and the four-operator expectation values 
including dissipation, we will in general have to find equations of motion of the type (Q). Taking into account the 
interaction with the laser, the equations for the different expectation values will be coupled. The task would then 
be to find the equilibrium expectation values, to which the observables will evolve at large times, as e.g. mP q in the 
case of eqn. (^). As we are here interested in the zero-temperature case only, we expect these equilibrium limits to 
be zero in general. It should then be justified to take the simpler route and use a regression theorem, i.e. the basic 
equations for the Fock operators, like (||), in conjunction with Wick's theorem, to calculate the time dependence 
of the expectation values. Such a procedure is readily justified in the case without decay. We will have to show, 
however, that such a regression theorem may also be used for the dissipative case, and give precise rules for its 
application. 

We start with the dynamic equation for the operators involved in the process to be studied. In the homogeneous 
limit, the quasiparticle modes are denoted by their momenta. We recall that the fluctuation field operator ( |l9| ) 
may be expanded in terms of particle Fock operators, -0(r, t) — V^ 1 ^ 2 ^2 k b k exp(ifc ■ r), and the particle operators 
may be expressed by quasiparticle ones through bk = UkPk + v k (3*_ k (cf. (fl9|)). The quasiparticle transformation 
coefficients are given by u k = (1 — a 2 ) -1 / 2 , v\ = u\ — 1, with a k = 1 + fc 2 — fc(2 + fc 2 ) 1 / 2 , fc — \k\/k Q , k = y/S-rrano, 
no being the condensate density. The Bogoliubov quasiparticle frequencies are Lu k — (hk 2 /2m) fc\/ 2 + fc 2 . 

We assume two-photon resonance, where the difference Ail = — fl s of the frequencies of the laser, fi;, and of 
the scattered photon, Q s , equals the frequency of the recoiling atom mode, Afi = u A k, with Afc = fc; — k s being 
the recoil momentum. If we neglect non-resonant terms, the coupled equations of motion read 

± ( P Ak {t) \ _( -lAk/2 -ih \ ( p Ak (t) \_( F Ak (t) \ 

dt \ at(t) J [ ztt ) \ dt(t) J { J> 

'0Ak(t) = -(7Afc/2)j9 A fc(t) - F_ Ak (t), (75) 

Here we have chosen an interaction picture by defining (3 q (t) = q (t) exp(iLJ q t), a s (t) = a s (t) exp(iCl s t), F q {t) = 
F q (t) exp(iu) q t), and fl = Cl(t) exp(if2;t). We denote by 17 = — 2|d ge • £ <+) (fci)| g ge (k,2)/A. the two photon Rabi 
frequency, given in terms of the dipole momentum d ge between the internal ground and excited states, of the 
positive frequency part of the electric field describing the laser mode, the coupling g ge (k2) between the scattered 
photon mode and the atom, and the detuning A from the excited state. Neglecting the decay width 7a& of 
the quasiparticle mode as well as the corresponding Langevin operator F A j., we may straightforwardly calculate 
the time evolution of the particle occupation numbers (n^ ), (ftj-Afc)> an( ^ °^ ^ ne squeezing parameter £3 = 
[(An°) 2 + (An^ fc ) 2 — 2 (( n®h h Ak ))]/( n" + n b Ak ), using their expressions in terms of quasiparticle operators given 
above: 

(n;> = (o a ot)-l, (76) 
(n b ±Ak ) = u 2 Ak ($l A J ±Ak )+v 2 Ak (($l Ak $ TAk ) + l) , (77) 
(A<) 2 = « >« <> + !), (78) 
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(An 



Ak) 



n b A k ', 



L Ak 



l Ak 



(0 

\(a s (3. 



'Ak , 



(79) 
(80) 



To be able to include collisional loss, however, a more careful approach is required. In appendix |A| we show that 
the dynamics of the expectation values appearing on the right hand sides of ([76 77 |8^) is given by the system of 
equations 



(a s al){t) + n° Ak 
(a s pAk)(t) - cc 



= exp 



I - 



lAk 

0^ 
2iVL -2ift - 



(a s (3 Ak )(t) -cc. / [\^-2ifi -2iQ -jAk/ 

^L A fc^-A*(*) > - Ak = exp[- 7 Afc(i - t )} ({pl A J- 




(t-t ) ( 




-n° 



PijAk)(to 
' a s al)(t ) - 
a s PAk){to) - cc 



Ak 



n Ak 



(81) 
(82) 



Note that the operator ordering in the expectation values of the density operators is crucial. This ordering may 
of course be changed in the final result using the commutation relations. Thus ( |s"l|) shows that applying the basic 
equations of motion (|7^ ) for the photon operators, at temperature T = 0, where the equilibrium quasiparticle 
population mP q vanishes, is only possible for their anti-normal-ordered products. More generally, products of the 
operators a s , f3 Ak , and their respective hermitian conjugates have to be normal ordered in the sense, that a s has 
to be interpreted as a creation operator and a\ as an annihilator, whereas the quasiparticle operators are taken 
in the usual sense. This is due to the fact that our basic Hamiltonian in the approximation which leads to the 
equations ( f74|) is not excitation number conserving due to its pair production term cx a\f3 Ak . For this reason we 
have already written in eqns. (|76|), (\n\j, and ( |S0|) the crucial operator products in the required ordering. 

The next step would be to derive corresponding equations for products of more operators, which can be done 
straightforwardly along the same lines. However, as we are for the first interested in the zero-temperature case, 
we can restrict the effort by following the operator ordering rule stated above and using the basic time evolution 
equations (f74|), (f75|). In appendix |a| we show that this rule, applied to products of an even number of operators, 
in the Markov approximation leads to the correct time evolution of the expectation values at T = 0. 
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FIG. 2. Time dependence of the mean particle occupation numbers of the scattered photon mode, (n£ 2 ) (dashed line), 
and of the modes of the recoiling atoms, ( n b Ak ) (solid line), { 'h b _ Ak ) (dashed-dotted line). The corresponding curves without 
decay are shown as dotted lines. 



In the rest of this paper we present the results for the time evolution of the occupation numbers and squeezing in 
the same system as considered in j^], but including quasiparticle damping. We consider specifically a homogeneous 
condensate of N = 10 6 23 Na atoms within a volume V = 10~ 8 cm 3 , which we choose such that the condensate 
density no = 10 14 cm~ 3 is of the order which is typically observed in experiments. With a — 2.8 nm we thus have 
ko = 2.65 • 10 6 m~ 1 , fik'^/2m — (2ir) 1.5 kHz. The time dependence of the occupation numbers and the squeezing 
parameter is shown in Figs. and [| These results have been discussed in detail in ||. We plot the logarithm of the 
occupation numbers ( ), ( h b Ak ), ( n b _ Ak ), for Ak = Afc/fco = 5, i.e. for coupling into the particle regime of the 
Bogoliubov spectrum, where oj k ~ frk 2 /2m. We have chosen a two photon Rabi-frequency of £1 = Is -1 . Initially 
the atomic occupation is (n b Ak (t = 0) ) = v Ak due to collisional interactions. Only for t > 0.1ms the photon and 
atom number occupations are of the same magnitude, (h k ) ~ ( n A/c)- For {fi h ± Ak ) approaching ./V = 10 6 , the 
calculation yields non-physical results since we have neglected the depletion of Nq. The range of times where the 
calculation is valid may however be increased by choosing a smaller Q. 
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The effects of the quasiparticle damping due to the finite width 7a/c become visible for t > 1ms. The curves 
shown correspond to a width of 75fc = 2.8 • 10~ 3 W5fc o = 7.0 • 10 2 s _1 . The damping affects all occupation numbers 
shown, also that of the photon mode, since the photon production rate depends on the occupation of the recoiling 
atomic mode. A hypothetical stronger decay rate would result in a time evolution, where ( n±/± k ) approach the 
initial value v Ak for large times and ( n k2 ) grows at a considerably lower rate. 




-6 -5 -4 -3 -2 -1 

log 1() (t/sec) 

FIG. 3. Time dependence of the squeezing parameters £i, i = 1, 2, 3, with (solid lines) and without quasiparticle decay 
(dashed and dash-dotted). 

Fig. ^ shows the squeezing parameter £3. For times t < 10 fis it is close to 1, i.e. there is essentially no squeezing 
compared to the classical case of independent modes. This is because the relative atom-photon number variance is 
still dominated by the initial "thermal" variance of the atom mode population, [An b Ak (t = 0)] 2 = ( n b Ak )(( h b Ak ) + 
1) = v Ak (v Ak + 1), which, in the particle regime where v Ak <§C 1, resembles the variance of a coherent state. (For 
recoil momenta in the quasiparticle regime, AA; < fco the initial value of £3 would be considerably larger than 1.) 

For t > 10 /is the squeezing starts to improve due to the atom-photon pair production, but only until t ~ 0.1 ms. 
It becomes clear that the quasiparticle damping causes the squeezing to deteriorate already at a considerably 
earlier time and larger value of £3 compared to the loss-less case. Decoherence due to quasiparticle damping is 
therefore a further important limiting factor. 

In the context of spin squeezing, we found analytically that within the restricted Fock space of the modes ±Afc, 
fc 2 , we have ( Ji ) = ( J 2 ) = for all times. Thus the mean abstract spin vector ( J ) points in 3-direction. Fig. || 
shows, that although the squeezing of all three spin components becomes worse for large times compared to the 
case of coherent modes, there is a relative squeezing of up to 2 orders of magnitude of the 3-component compared 
to the 1- and 2-components which are of equal size 

V. CONCLUSIONS 

In this paper we have studied the dynamics of interacting condensed Bose-gases at finite temperature. We have 
set up a system of equations of motion for the operators and operator products which describe the fluctuations 
around the mean field and used the linear response approximation to describe the time evolution of the system. We 
have focused on the dissipation in the time evolution of systems which initially are out of equilibrium, specifically 
for initial states with a single excited (quasiparticle) eigenmode. As is well known, the decay of this mode is 
caused predominantly by interaction processes of the quasiparticles which involve a condensate atom and which 
can be classified as Landau or Beliaev processes. We have reviewed the derivation of the general expressions of the 
quasiparticle decay rates for a trapped condensate at finite temperature and extended this to the case where the 
atoms in the initially excited mode are in a different internal state than the condensate atoms. The interactions 
taken into account in this two-level case are collisions between atoms in different internal states which conserve the 
internal quantum number. We calculated the coupling functions which enter the expressions of the Beliaev and 
Landau decay rates and found the same infrared ^-proportionality of the Beliaev rate as for the single internal 
level case. 

In order to gain a consistent description of the dissipative dynamics of densities and higher correlation functions 
in the long time limit we used the Heisenberg-Langevin approach. In this way the time evolution equations of 
correlation functions including their correct equilibrium limit may be determined straightforwardly. In this paper 
we are interested in the zero temperature dynamics of expectation values, which yield the number squeezing 
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between different Fock modes of a system. We were able to introduce and justify a regression theorem which we 
used to compute higher correlation functions using the basic equations of motion for the Fock operators only. 

Specifically, we considered Bragg-scattering of photons off a finite homogeneous Bose-Einstein condensate of 
sodium atoms and studied the relative occupation number squeezing between the scattered photon and the recoiling 
atomic mode. The parameter which quantifies the squeezing measures the variance of the number difference of the 
modes compared to its classical limit which is given by the sum of the mean occupations. Hence our calculations 
required to determine the evolution of densities as well as mean values of products of four operators. Our numerical 
results show that there is a substantial limitation of squeezing due to the collisional interactions between the atoms. 
Compared to the approximation where quasiparticle interactions have been neglected the achievable squeezing is, 
in addition, considerably reduced by decoherence due to dissipation. 
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APPENDIX A: OPERATOR ORDERING RULE 

In this appendix we derive the time evolution equations (jsi|), (|82]), for the expectation values appearing on the 
right hand sides of ( f76||77|) . Moreover we justify our procedure, that, for T = and in the Markov approximation 
we may compute expectation values of normal ordered products of more than two operators by using the basic 
equations ( |74| , |75| ), without taking into account the Langevin operators. Normal ordering here means, that a s is 
treated as a creation operator, and a\ as an annihilator, whereas the quasiparticle operators are treated in the 
usual way. 

In order to obtain the dynamical equations for the expectation values of the required products of two operators, 
we need to solve the following equations which are derived from ( |74] , |75| ). 



d_ 

db 



(a s 5,t) = 0^ 0^ in (a s a\) - ). (Al) 

i a s j3 Ak ) - c.c. / \-2ift -2iti - lAk /2 J \ { aJ Ak } - c.c. / \ (F Ak al\ 

J t (pl A J-Ak) = -lAk(Pl A J-Ak)- ({FtJ-Ak) +C.C.) . (A2) 

The expectation values involving the Langevin operator determine the equilibrium values of the observables. Using 
the formal integration of (|74|), together with eqn. (|^) for the two time correlation function, as well as eqns. (|2§), 
(p3|), (Ell) (which show that the expectation values of the Langevin operator with the Fock operators at time to 



vanish), we find: ( P]_ Ak (3± A k ) = TAfe^Afc-i ^Akfit ) = ^ This leads, upon integration of (Al), (|A|), to the final 
result for the time evolution, eqns. (|8l|), (p2|). 

The crucial condition for the application of the regression theorem is the operator ordering in the expectation 
values of the density operators. The general rule is that products of the operators a s , /?Afej and their respective 
hermitian conjugates have to be normal ordered. This is meant in the sense, that a s has to be interpreted as a 
creation operator and a\ as an annihilator, whereas the quasiparticle operators are taken in the usual sense. In 
the following we sketch the proof that using the basic time evolution equations ( |74| , |75| ) and Wick's theorem for 
products of an even number of operators, in the Markov approximation leads to the correct time evolution of the 
expectation values at T = 0. 

The time evolution equation of a (in the above stated sense) normal ordered product of n Fock operators, using 
the chain rule and eqns. (74 7q), involves terms containing the Langevin operator. It is the expectation values of 



these terms which we have to consider and which need to vanish for the desired procedure to be justified. Upon 
insertion of the formal solutions of eqns. ([M 75) for the Fock operators in these terms, their value is determined by 



expectation values of mixed products of Langevin operators at different times and Fock operators at time to ■ The 
normal ordering of the initially considered operators carries through to a normal ordering of these mixed products. 
Moreover, the terms involve time integrals over those Langevin operators which entered through the integrated 
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eqns. (f74j), (|7^)- For the case of four operators which we need to compute squeezing, these integrals are e.g. of the 
form 

f dt'dt"dt"'(F\t)F\t')F(t")F(t"')), 

[ dt'(F\t)a\t )a(t )F(t')), (A3) 

where a is either of the operators a\, /3±Afc- When inserting ( [4l| ) we see that, at T = 0, only those terms can be 
different from zero, which contain an operator gij(to) at the very left and one gtAto) a t the very right. (If there 
is an a s (to) at the very left and/or an aj(io) at the very right, the expectation value may be broken down using 
Wick's theorem, and the problem reduced to the case of n — 2 or less operators.) 

The final step is to show that the remaining non-vanishing terms are zero in the Markov approximation. This 



is done as before by assuming that the coupling functions (pq-|2q) which multiply the integrals (A2) vary only 
slowly over the range of modes, such that the sums over the mode indices result in functions which are strongly 
peaked at a single time and may be written as proportional to delta-distributions in time. Integration over the time 
argument of one of the Langevin operators at the very right or left then immeditately yields, that the remaining 
terms are all proportional to S(uJAk + + u>j). Since we consider | Afc| ^ only, the argument of the distribution 
is positive. The expectation values in the equations of motion, which involve Langevin operators, therefore vanish 
in the Markov approximation at T — 0, for the same reason as this is the case, with n° q — 0, for (§^). 
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